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Abstract. 

N-body simulations of collisionless collapse have offered important clues to the construction of realistic stellar dynamical 
models of elliptical galaxies. Understanding this idealized and relatively simple process, by which stellar systems can reach 
partially relaxed equilibrium configurations (characterized by isotropic central regions and radially anisotropic envelopes), is a 
prerequisite to more ambitious attempts at constructing physically justified models of elliptical galaxies in which the problem 
of galaxy formation is set in the generally accepted cosmological context of hierarchical clustering. 

In a previous paper, we have discussed the dynamical properties of a family of models of partially relaxed stellar systems (the 
/ w models), designed to incorporate the qualitative properties of the products of collisionless collapse at small and at large radii. 
Here we revisit the problem of incomplete violent relaxation, by making a direct comparison between the detailed properties 
of such family of models and those of the products of collisionless collapse found in A'-body simulations that we have run for 
the purpose. Surprisingly, the models thus identified are able to match the simulated density distributions over nine orders of 
magnitude and also to provide an excellent fit to the anisotropy profiles and a good representation of the overall structure in 
phase space. The end-products of the simulations and the best-fitting models turn out to be characterized by a level of pressure 
anisotropy close to the threshold for the onset of the radial-orbit instability. The conservation of Q, a third quantity that is argued 
to be approximately conserved in addition to total energy and total number of particles as a basis for the construction of the f (v) 
family, is discussed and tested numerically. 

Key words, stellar dynamics — galaxies: evolution — galaxies: formation — galaxies: kinematics and dynamics — galaxies: 
structure 



1. Introduction 

The collapse of a dynamically cold cloud of stars can lead to 
the formation of realistic stellar systems, with pro jected density 
profiles well represented by the R 1 ^ 4 law (van Albada 1982). 
The theoretical framework for the mechanism of incomplete 
violent relaxation th at governs this proce ss of structure forma- 
tion was proposed by Lvnden-Bell ( 1967), who argued that fast 
fluctuations of the potential during collapse would lead to the 
formation of a well-relaxed isotropic core, embedded in a ra- 
dially anisotropic, partially relaxed halo. This general picture 
served as a physical justification for the construction of the so- 
called /oo models, which indeed recovered the R 1 ^ 4 law and, 
suitably extended to the case of two-component systems (to 
account for the coexistence of luminous and dark matter), led 
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to a number of inte resting applications to the observations (see 
iBertin & Stiavellil 1 984l fl993l and references therein). 

An attempt at deriving the relevant distribution function 
directly from the statistical mechanics of incomplete vio- 
lent relaxation suggested that, in addition to the /«, models, 
one could consider alternat ive models, called the / (v) mod- 
els (Stiavelli & Bertin 1987), with similar overall characteris- 
tics. The key ingredient for the construction of the /^ v) dis- 
tribution function is the conjecture that a third quantity Q, in 
addition to the total mass M and the total energy tLj ^ is ap- 
proximately conserved during the process of collisionless col- 
lapse (of course, we are referring to systems characterized by 
vanishing total angular momentum, J tot = 0). This quantity 
is introduced to model the process of incomplete violent re- 
laxation, ensuring a radially biased pressure tensor and a 1/r 4 
density profile in the outer parts of the system. Because of 
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their relatively straightforward derivation from the Boltzmann 
entrop y, these models were revisited recently (iBertin & T renti 
2003) and used to demonstrate the onset o f the gravother- 
mal catastrophe dLvnden-Bell & Woodlll968l) for such a one- 
parameter sequence (at fixed v) of anisotropic equilibria; a pre- 
liminary inspection of the general characteristics of the / (v) 
models then convinced us that, with significant advantage over 
the /oo models, they might also serve as a good framework to 
interpret the results of simulations of collisionless collapse not 
only qualitatively, but also in quantitative detail. Therefore, we 
proceeded to examine their intrinsic properties systematically 
jTrenti & Bertinl200 4: hereafter Paper I), and we will take ad- 
vantage of that work for the study presented in this paper. 

In this article we describe the results of a relatively wide 
set of numerical simulations of collisionless collapse, aimed 
at studying the phase space evolution and settling of the sys- 
tem during violent relaxation, and we then compare in detail 
the properties of the quasi-equilibrium end-products thus ob- 
tained with those of the / (v) models. In particular, we discuss 
the role played by the initial conditions and find that a certain 
degree of dumpiness is required for an efficient mixing in the 
single-particle angular momentum distribution; this form of re- 
laxation turns out to be crucial for a good match with the 
family of models. The Q conservation is then studied directly 
by looking at its time evolution. For a significant range of col- 
lapse factors, as determined by the initial values of the virial 
ratio u - {2K/\W\), = q, an approximate conservation is indeed 
observed. The end-products (and thus the best-fitting models) 
tend to be characterized by a value of the global anisotropy 
parameter close to marginal stab ility with r espect to the radial 
orbit instability ( Polvac henko & Shukhmanll98ll) . 

The paper is organized as follows. After introducing our ba- 
sic models and notation (Sect. [5), in Sect.[3]we start by review- 
ing the choice of the numerical code and we then chec k the re- 
sults obtained in some test runs against the tree code o flDehner] 
(200(3). In Sect.|4]we discuss the initial conditions adopted for 
the simulations of collisionless collapse, with special attention 
to the issue of dumpiness in phase space. In Sect. [5] we charac- 
terize the end-products of the simulations in terms of a few key 
indicators (i.e., central concentration, global anisotropy, den- 
sity and anisotropy profiles, deviations from spherical symme- 
try) and describe their dependence on the initial conditions. In 
Sect. |6]we examine the hypothesis of the approximate conser- 
vation of Q. We then move, in Sect. to the comparison of 
the end-products of the simulations with the / (v) models (in 
terms of density and anisotropy profiles and directly in phase 
space). In Sect. [8] we draw the main conclusions from this 
study. Finally, in the Appendix we provide additional com- 
ments on the issue of dumpiness in phase space. 

2. f (y) models, units, and notation 

In general, we will keep the same notation as in Paper I. We 
recall that the relevant distribution function is obtained by 
extremizing the Boltzmann entropy S - -Jflogfd 3 xd 3 w 
at fixed total mass M - f fd 3 xd 3 w, total energy E lo , - 
1/3 f Ed 3 xd 3 w, and Q = J J v \E\- 3v/4 fd 3 xd 3 w. Here E and 
J denote single-star specific energy and angular momentum, 



while x and w denote positions and velocities respectively. 
This leads to the function f v) = A exp [-aE - d(J 2 /\E\ 3/2 ) v/2 ], 
where a, A, d, and v are positive constants. The / (v) models are 
then constructed by solving the Poisson equation for the un- 
known potential <t>(r) numerically. At fixed value of v, one may 
think of the free constants as providing two dimensional scales 
(for example, M and E tot ) and one dimensionless parameter, 
such as = -aO(O), the central depth of the dimensionless 
potential well. By (1; 5) / (v) model we will denote the model of 
the f {v) family with v = 1 and *P = 5. 

The / (v) models represent equilibrium configurations de- 
signed to describe the products of incomplete violent relax- 
ation. They are characterized by a density profile p(r) falling off 
as 1 /r 4 at large radii and as 1/r 2 in the inner part of the system, 
outside a central "core". The size of the core becomes smaller 
as the concentration parameter ^ increases. On the large scale, 
apart from such freedom in central concentration and core size, 
the shape of the density profile is basically independent of the 
(v; W) parameters (see Fig. 3 in Paper I). Interestingly, although 
this feature had not been imposed at the beginning (when the 
function / (v) is constructed), the projected density distribution 
of the / (v) models is typically well fitted, on the large scale, 
by the R 1 ^ 4 law; residuals in the fit are reduced if one consider s 
the generalized R 1 ^" law (with n a free parameten lSersid f9 68 ). 
depending on *P (see Figs. 4-5 in Paper I). 

In contrast with other approaches (e.g., see Osin kolfl979l 
and lMerrit3l 1 985l) where the anisotropy profile is assigned a 
priori, in the / M models the velocity dispersion anisotropy pro- 
file a(r), defined as a(r) = 2 — ((wjj) + (w^))/(w 2 ), must be 
computed a posteriori and its shape depends on v and *P (see 
Fig. 6 in Paper I). The structure of the distribution function only 
guarantees that the models match the asymptotic requirements 
suggested by the picture of incomplete violent relaxation, i.e. 
at large radii, where the pressure is radial, and in the central 
regions, where the pressure is isotropic. The global anisotropy, 
measured by the quantity 2K r jKj, i.e. twice the ratio of the ra- 
dial to the tangential kinetic energy, depends on the choice of 
(v;*F) and correlates with the central concentration (e.g., see 
Fig. 7 in Paper I). Models with *P < 4 are characterized by an 
excessive degree of radial anisotropy (i.e. 2K r /Kj > 1.7), and 
are thus unstable. 

The physical system of units adopted in this paper is de- 
fined by 10 kpc for length, 10 11 M Q for mass, and 10 s yr for 
time. In this system, natural for studies on the galactic scales, 
velocities are measured in units of « 97.8 km/s and the value 
of the gravitational constant G is 4.497 1 . 

The majority of simulations consists of runs starting from 
20 cold clumps of 16 kpc radius in a sphere of 40 kpc radius, 
with u = (2K/\W\) t =o in the range 0.05-0.25. After the col- 
lapse the system has a half-mass radius around 8 kpc. The to- 
tal mass of the system is 2 x 10 11 M G . The dynamical time, 
which we define as tj = GM 5 ' 2 /(—2E tot ) 3 ^ 2 , is therefore typi- 
cally » 1.2 x 10 s yr, i.e. 1.2 in our units. As a result, when we 
stop the simulation at time 80, the system has evolved for sev- 
eral tens of dynamical times. In any case, we should recall that 
the results obtained are scale-free, that is they can be rescaled 
to other choices of mass and radius if so desired. 
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3. The code 

Direct N-body simulations of self-gravitating stellar systems 
require huge amounts of computing time because of the TV 2 
scaling of the code complexity with the number of particles 
employed. To model the evolution of collisionless systems, 
several algorithms have been developed that treat the gravi- 
tational interactions approximately, with a lower complexity. 
In principle, to study the process of collisionless collaps e we 
have two options: either a tree code (e.g., [Barnes & Hut 1986) 
or a particle-mesh like algorithm (e.g., see Ivan Albadal ll982: 
iMcGlynnligsllHernauist & OstrikeJl992h . Tree codes are in- 
trinsically slower, scaling as N log N, with respect to particle- 
mesh Poisson solvers, for which the complexity is linear in N, 
but the latter have the disadvantage of a lower spatial resolu- 
tion, being limited by the size of the grid used. 

In this paper we are interested in the large scale structure of 
the end-products of collisionless collapse, for systems that do 
not exhibit large deviations from spherical symmetry. The natu- 
ral choice thus appears to be that of a particle-mesh code, based 
on a spherical grid and an expansion in spherical harmonics. 

The code used in t he present study is thus a new version of 
the Ivan Albadal Jl982n code. The relevant changes introduced 
are briefly described below, in Sect. 13. II For completeness, we 
have also run (see Sect. 13. 21 a numb er of compariso n simula- 
tions with the fast code developed by DehnerJ J2000l) . 

3.1. An improved particle-mesh code 

The key feature of the Ivan Albadal dl982h code is the solution 
of the Poisson equation V 2 <I> = AnGp, which relates the mean 
potential cD of the system to the mass density p, by means of 
Fourier techniques. Once the potential has been computed by 
expanding the density in spherical harmonics, the acceleration 
is obtained by numerical differentiation, and the particles are 
advanced by a fixed time step, using a leap-frog scheme. 

At variance with the original implementation, to preserve 
accuracy and to avoid systematic errors, we decided to drop the 
angular grid and to treat the angular dependence of the force ex- 
actly, in terms of the single-particle Legendre polynomials (for 
further details, see lTrentl2004 . Basic ally, this choice c hanges 
the code in the direction of the code of McGlvnn (1984) and of 
the self-consistent field code of Hernquist & Ostriker ( 1992). 
We preferred to keep the radial grid because of its flexibility 
with respect to the density profile, especially under conditions 
of rapid evolution (we use a subroutine to generate a grid con- 
taining a fixed fraction of the total mass in each shell). 

The density is assigned to the radial grid by means of a 
cloud-in-cell scheme with a linear kernel, i.e. a particle con- 
tributes to the density of the two closest cells with a weight 
depending linearly on the distance from the center of the cell 
considered. The same kernel is then used to assign the force 
from the grid to the particle. The time step is adaptively chosen 
in such a way that particles are not allowed to cross more than 
one radial cell during one step. 

The code has been tested extensively, in terms of its accu- 
racy in conserving total energy and total angular momentum 
for equilibrium and non-equilibrium initial conditions and in 
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Fig. 1. Total energy conservation (left) and angular momentum 
(right) for a simulated ( 1 ; 5) equilibrium model with 4 x 1 5 
particles. The mass of the system is M — 1, the total energy 
Ewt = -1.73, and the half-mass radius tm = 0.5. The magni- 
tude of the total angular momentum associated with the initial 
conditions corresponds to a value of the dimensionless rotation 
parameter A = J,JE tot \ l/2 /(GM 5/2 ) * 1(T 4 . We recall that time 
is given in units of 10 8 years, which, in this case, corresponds 
to * 1.4 t d , with t d = GM 5/2 /(-2E tot ) 3/2 . 

conserving single-particle energy and angular momentum for 
runs of spherically symmetric equilibrium models. The mod- 
ified Poisson solver scheme combined with the adaptive grid 
ensures a significantly better accuracy with respect to the orig- 
inal implementation. The typical total energy and total angular 
momentum conservation for a run with 10 5 particles is of the 
order of 10~ 5 per dynamical time in quasi-stationary configu- 
rations (see Fig.[Q. 

3.2. Comparison with Dehnen's code 

As a further test, we have also run a few test simulations by 
comparing, under identical conditions, the perf ormance of our 
code to that of the fast tree code Gvr FalcON jDehnenl F2000 . 
2002), within the NEMO distribution llTeubenlll995h . In such 
tests, we adopt the following procedure. We first generate the 
initial conditions in the physical units used by our code (see 
Sect. |2J and we run the simulation. We the n convert the initial 
cond itions to the natural units defined by iHeggie & Mathieul 
(1986) and used in Dehnen's code. Finally, we run the simu- 
lation within the NEMO environment. The quality of the inte- 
gration is checked with the standard NEMO tools of analysis. 
At the end of the simulation, a "snapshot" of the system is ex- 
ported and converted back to our units, in such a way that it 
can be processed by the same diagnostics used for the particle- 
mesh code. 

The initial conditions for these runs have been chosen in 
order to be representative of the sample investigated in this pa- 
per; they are described in Table[0(C4.1 and C4.3 entries), with 
the properties of the final equilibrium state listed in Tabled 

For the runs with Dehnen's tree code we adopted the fol- 
lowing choice of integration parameters: tolerance parameter 
9 = 0.5 (standard choice 0.6) to improve accuracy in the cal- 
culation of forces; softening length e = 0.01 (in natural units; 
standard choice 0.05) to increase central resolution; minimum 
allowed time step 1/2 8 (i.e. « 724 steps per dynamical time). 
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Fig. 2. Com parison between our code and GyrFalcON 
llDehnenl2 000) for the run C4.1. Final density (upper left) and 
anisotropy (bottom left) profiles, and single-particle energy dis- 
tribution (bottom right) for a run starting from 10 5 particles in 
10 cold clumps withM = 0.15 (see Tables[2and|2]f or further de- 
tails about the simulation); here the solid lines give the results 
obtained with GyrFalcON, while the crosses refer to our code. 
In the upper right panel we plot the differences in the density 
profile Ap = 2(p fakon -Ppm) lip falcon +P P m) (crosses) and in the 
anisotropy profile Aa = af a k-on - apm (solid line). 

With this choice of integration parameters, the energy and an- 
gular momentum conservation is very good: in one dynami- 
cal time tj, the relative changes are AE tot /E tot < 10~ 5 and 
AJ tot /J tot < 10- 4 . 

The required CPU time to complete the simulation is 
marginally higher than with our code, which, however, has not 
yet been optimized for speed. 

As desired, for these runs we find a substantial similarity 
in the properties of the end-products obtained by the two dif- 
ferent methods of integration (see Fig. |2j. To be sure, small 
differences naturally arise, as expected. The main systematic 
difference is in the degree of anisotropy characterizing the end- 
products of the simulations. In fact, the output from the tree 
code is slightly more isotropic: the final global anisotropy, mea- 
sured by 2K r /K T , is up to 7% lower, with a slight outward shift 
of the anisotropy profile, corresponding to a more efficient core 
relaxation. This might be related to the residual collisionality 
present in the tree code. 

4. Choice of initial conditions 

If the initial conditions are not too artificial, during the pro- 
cess of collisionless collapse, violent relaxation can take place, 
with significant mixing in phase space, and wipe out much 
of the details that characterize the initial conditions. In real- 



ity, violent relaxation is incomplete. Therefore, the final state 
is that of an approximate dynamical equilibrium character- 
ized by an anisotropic distribution function, different from 
a Maxwellian (which would correspond to thermodynamic 
equilibrium). Because of such incomplete relaxation, the end- 
products of the simulations do conserve some memory of the 
initial state. 

4.1. Uniform initial conditions, clumpy initial conditions, 
and the cosmological framework 

Some of the papers addressing the problem of collisionless col- 
lapse start from "uniform" initial conditi ons in pos ition and 
velocity space. For example Aguil ar & Merrittl |l990) assume 
an initial 1/r density profile and then explore the way the col- 
lapse proceeds, by varying, in addition to the initial virial ratio 
u = (2K/\W\), = o, the shape of the initial density profile (by 
shrinking t he syst em along one axis) and the amount of rota- 
tion. llJdrvl Jl993l) starts from uniform cold spheres, and also 
varies, in addition to the above-mentioned parameters, the ini - 
tial anisotropy content 2K r /Kj. Recently, Boil v et al 
starting from cold uniform spheres or spheroids, focus on the 
effects introduced by the number of particles used in the simu- 
lation. 

A few earlier investigations (I van Albadalll982t iMcGlvnnl 

Il984t iMav & van Albadal Il984l lLondrillo et all Il99ll) cori> 
pared "clumpy" to "uniform" (or "homogeneous") initial con- 
ditions, showing that clumpy initial conditions lead to end- 
states with pr ojected density d istributions well fitted by the R 1 ^ 4 
law (although Aguilar & Merritt 1990 point out that, for very 
small values of u, the R 1 ^ 4 law is approximat ely rec overed even 
for homogeneous initial conditions). I Udrvl dl993h argues that 
starting from a multi-component initial mass spectrum for the 
simulation particle distribution can be be an alternative way to 
represent a clumpy initial density configuration. However, the 
introduction of simulation particles as massive as to be repre- 
sentative of clumps, would introduce effects of dynamical fric- 
tion that per se would go beyond the picture of collisionless 
violent relaxation.] 

Recently IRov & Perezl J2004I) studied the outcome of vio- 
lent collapse starting from an initial uniform background with 
the possible addition of small clumps of stars. Although their 
clumpy initial conditions are rather different from those con- 
sidered in the present paper, they also noted that clumpy simu- 
lations lead to steep density profiles with small cores. 

As will also be demonstrated later on (see Sect. l5.3> . the key 
point that distinguishes clumpy from uniform initial conditions 
is that, in general, only the former allow for significant mixing 
in phase space, thus making it possible for violent relaxation to 
proceed properly. Thus in this paper we will focus on simula- 
tions starting from clumpy configurations. As discussed in the 
next Subsection, for the present study we do not require that 
our initial clumps be in internal dynamical equilibrium, since 
their purpose is to avoid excessive homogeneity in the (E, J 2 ) 
phase space (see also Appendix). In particular, the clumps are 
not intended to be a realistic representation of possible condi- 
tions at a given epoch in the past. In fact, the effects of violent 
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relaxation become important in a few dynamical times, inde- 
pendently of the precise epoch when the process is imagined to 
occur. 

To be sure, to identify a realistic set of initial conditions, 
one should consider a satisfactory cosmological framework. 
We plan to do this in future investigations, because this would 
lead us well beyond the scope of the present paper. In this 
respect, the use of clumps is already one important step for- 
ward with respect to the use of homogeneous initial conditions. 
Eventually, we should devise a method for determining a "spec- 
trum" of clumps with properties compatible with the expecta- 
tions of current cosmological scenarios (see also lKadfl99ll 
and further discussion in the Appendix). For the moment, we 
are satisfied with identifying the initial conditions under which 
sufficient mixing in phase space is guaranteed. Note that cos- 
mologically oriented simulations are centered on the clustering 
and growth of dark matter halos, while in this paper, given our 
focus on the R l/4 law and on the deviations from it, we have in 
mind mostly luminous matter. 

4.2. Setting up clumpy initial conditions 

In a clumpy initial state the N particles are grouped in Nc 
spherical clumps, each of them containing Nj stars, so that 
N = 2^ Nt, with (Ni) = N/Nc- Within each clump the star 
distribution is homogeneous. The centers of mass of the clumps 
are distributed uniformly inside a sphere of radius R, which 
defines the size of the system at the beginning of the simu- 
lation. The clump radius is Rc, with Rc < R, but such that 
Nc x R c > R 3 (this condition ensures that the sphere of ra- 
dius R is well filled by stars). The initial kinetic energy may 
be associated with the ordered motion of the center of mass of 
each clump (this is our default choice for the simulations of 
type C described below; in this case the velocity is assigned 
by drawing from an isotropic distribution) or with the random 
motions of the stars within the clump (in these cases we add a 
subscript h to the simulation label; here the center of mass of 
each clump is taken to be at rest). In general, with this choice 
of initial conditions the clumps are not in internal dynamical 
equilibrium. We note that when the number of clumps used 
is low, the initial configuration may deviate significantly from 
spherical symmetry (with projected shapes up to those of an E3 
galaxy). Formally the limits Nc — > N and Nc — > 1 both lead 
to the case of homogeneous initial conditions. 

In the case of homogeneous initial conditions (simulations 
of type U and S ), which we run for comparison, we employed 
two kinds of distribution: (1) a constant density within a sphere 
of radius R; (2) a symmetrized version of a given clumpy con- 
figuration (simulations of type S). The symmetrization process 
in (2) is performed by accepting the radius and the magnitude 
of the velocity of each particle, following the procedure for ini- 
tial clumpy conditions, but by redistributing uniformly the an- 
gular variables. 

In principle, we have a wide parameter space to explore, 
because we have to deal with the initial virial ratio u, the num- 
ber and size of the clumps, the cold/hot choice for the initial 
kinetic energy distribution (and the intermediate range of pos- 
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Fig. 3. Typical projected distributions in position (left) and ve- 
locity (right) space for hot (upper panels, run C4.1/,) and cold 
(lower panels, run C4.1) clumpy initial conditions. 



sibilities), the spatial distribution of the centers of mass of the 
clumps a nd of the stars within each clump. As noted earlier 
(e.g., see lLondrillo et aflll99ll) . we anticipate that the main 
controlling physical parameter is the initial virial ratio. 

Table^lists for each simulation the following information: 
the number N of particles used, the number of clumps Nc, the 
initial virial ratio u, the initial values of the shape parameters 
eo, rjo (based on the length of the axes of the homogeneous 
ellipsoid associated with the inertia tensor, taken to be in the 
order a > b > c, so that rj = c/a and e = b/a; the inertia ten- 
sor is referred to the particles within a sphere of radius 3tm), 
and the initial concentration C P0 - (p(0)/p(rM)), = o- As a sum- 
mary for the notation used, we note the following. We have 
divided the set of clumpy simulations in five subsets, from CI 
to C5. The simulations belonging to CI start with 10 5 particles 
in 10 clumps, the positions of which are fixed. The C2 series 
is a high resolution (8 x 10 5 particles) version of CI, but uses 
instead 20 clumps. In the C3 (high resolution, 8 x 10 5 parti- 
cles) and C4 (10 5 particles) series we use different seeds for 
the initial positions of the clumps and we also change other pa- 
rameters as described in Tabled Runs CV5.1 and CP5.2* are 
test runs especially performed to clarify some issues related to 
dumpiness (see Appendix). CV5.1 has clumpy conditions in 
velocity space as in run C4.1, but uniform homogeneous con- 
ditions in position space; in turn, CP5.2* has a clumpy con- 
figuration in position space (40 clumps of 6 kpc, with a filling 
factor Nc x R 3 C /R 3 = 0.135) and uniform conditions in veloc- 
ity space. Runs U refer to uniform homogeneous spheres (here 
the seed for the random numbers is not relevant given the high 
symmetry of the configuration) and the S series refers to the 
symmetrized runs. 
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Table 1. Initial conditions for the simulations. After the simu- 
lation identifier the columns list the number of particles N, the 
number of clumps Nc, the initial virial ratio u, the initial values 
of the shape parameters eo and 770, and the initial concentra- 
tion C p (). For the exact definitions and for the general charac- 
teristics of the groups (CI to C4, U and S), see description in 
Sect. 4.2. Four simulations, marked by a *, have been carried 
out with GyrFalcOn jDehnenll200fj|) : two of them (C4.1* and 
C4.3*) start from the same initial conditions as for C4.1 and 
C4.3 respectively, while C4.4* and CP5.2* are characterized 
by small clumps (with radius 2.8 kpc and 6 kpc respectively, 
distributed in a sphere of radius 40 kpc). 
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10 


0.1 


0.83 


0.70 


3.0 


C1.9 


10 5 


10 


0.075 


0.83 


0.70 


3.0 


CI. 10 


10 5 


10 


0.05 


0.83 


0.70 


3.0 


C2.1 


8 ■ 10 5 


20 


0.23 


0.93 


0.73 


2.8 


C2.2 


8 • 10 5 


20 


0.17 


0.93 


0.73 


2.8 


C2.3 


8- 10 5 


20 


0.12 


0.93 


0.73 


2.8 


C2.4 


8 • 10 s 


20 


0.06 


0.93 


0.73 


2.8 


C3.1 


8- 10 5 


20 


0.08 


0.95 


0.91 


2.2 


C3.2 


8- 10 5 


20 


0.18 


0.86 


0.80 


2.6 


C3.3 


8 ■ 10 5 


20 


0.15 


0.84 


0.70 


3.1 


C3.4 


8 • 10 5 


20 


0.23 


0.88 


0.73 


2.0 


C3.5 


8- 10 5 


20 


0.15 


0.95 


0.88 


3.7 


C3.6 


8- 10 5 


10 


0.15 


0.86 


0.80 


3.7 


C4.1 


10 s 


10 


0.15 


0.87 


0.80 


2.8 


C4.1* 


10 5 


10 


0.15 


0.87 


0.80 


2.8 


C4.1„ 


10 s 


10 


0.15 


0.87 


0.80 


2.8 


C4.2 


10 s 


20 


0.25 


0.75 


0.63 


1.5 


C4.3 


to 5 


80 


0.14 


0.90 


0.77 


1.9 


C4.3* 


10 s 


80 


0.14 


0.90 


0.77 


1.9 


C4.4* 


10 s 


80 + 


0.15 


0.85 


0.78 


2.0 


C4.5 


10 5 


400 


0.23 


0.99 


0.95 


0.8 


C4.5, 


10 s 


400 


0.23 


0.99 


0.95 


0.8 


CV5.1 


10 5 


10 


0.23 


1.00 


1.00 


1.0 


CP5.2* 


10 5 


40 


0.15 


0.81 


0.78 


0.7 


U6A 


8- 10 5 


N/A 


0.10 


1.00 


1.00 


1.0 


U6.2 


8- 10 5 


N/A 


0.19 


1.00 


1.00 


1.0 


U63 


8- 10 5 


N/A 


0.29 


1.00 


1.00 


1.0 


U6A 


8- 10 5 


N/A 


0.39 


1.00 


1.00 


1.0 


54.2 


10 s 


N/A 


0.25 


1.00 


0.99 


1.5 


54.3 


10 5 


N/A 


0.15 


1.00 


1.00 


1.9 



Table 2. Final configurations for the simulations of collision- 
less collapse listed in Tabled The column entries are described 
in Sect. |3 Note that the anisotropy profile in homogeneous sim- 
ulations can be non-mono tonic; this is indicated by 1 . All simu- 
lations of type CI and C2 start from identical initial conditions 
within each series, except for a constant scaling of velocities. 
The quantity AQ is referred to v = 5/8 in simulation CV5.1, to 
v = 1 in simulation U6A and to v = 3/4 in simulation U6.2; 
this is indicated by # . 





AM 


A2 


c P 


K 




e 


n 


CIA 


0.00 


0.13 


570 


1.61 


1.02 


0.91 


0.73 


C1.2 


0.002 


0.17 


600 


1.60 


0.94 


0.91 


0.74 


C1.3 


0.01 


0.20 


680 


1.59 


0.94 


0.90 


0.76 


C1.4 


0.01 


0.24 


790 


1.57 


0.88 


0.95 


0.79 


C1.5 


0.02 


0.30 


720 


1.52 


0.88 


0.96 


0.81 


C1.6 


0.03 


0.38 


820 


1.50 


0.93 


0.99 


0.80 


C1.7 


0.04 


0.44 


760 


1.47 


0.92 


0.97 


0.78 


CI. 8 


0.05 


0.52 


850 


1.53 


0.87 


0.96 


0.79 


C1.9 


0.06 


0.66 


1130 


1.67 


0.75 


0.97 


0.75 


CI. 10 


0.08 


0.72 


1090 


1.74 


0.79 


0.94 


0.69 


C2.1 


0.01 


0.13 


110 


1.52 


1.49 


0.87 


0.78 


C2.2 


0.02 


0.25 


160 


1.62 


1.24 


0.88 


0.78 


C2.3 


0.03 


0.4 


270 


1.70 


0.83 


0.81 


0.69 


C2.4 


0.07 


0.5 


520 


1.76 


0.74 


0.81 


0.63 


C3.1 


0.003 


0.47 


1690 


1.99 


0.44 


0.90 


0.73 


C3.2 


0.001 


0.26 


1250 


1.85 


0.55 


0.93 


0.70 


C3.3 


0.04 


0.57 


430 


1.60 


1.34 


0.92 


0.71 


C3.4 


0.02 


0.23 


500 


1.65 


1.15 


0.93 


0.81 


C3.5 


0.005 


0.24 


950 


1.73 


0.57 


0.96 


0.72 


C3.6 


0.005 


0.27 


690 


1.79 


0.75 


0.80 


0.77 


C4.1 


0.005 


0.27 


440 


1.77 


0.83 


0.80 


0.73 


C4.1* 


0.01 


0.27 


360 


1.68 


0.97 


0.89 


0.74 


CAA h 


0.00 


0.18 


240 


1.86 


0.51 


0.86 


0.83 


CM 


0.12 


0.10 


160 


1.40 


1.65 


0.90 


0.78 


C4.3 


0.10 


<0.01 


70 


1.60 


1.53 


0.84 


0.74 


C4.3* 


0.07 


0.15 


70 


1.50 


1.56 


0.86 


077 


C4.4* 


0.04 


0.4 


4000 


1.15 


5.30 


0.98 


0.96 


C4.5 


0.125 


0.02 


20 


1.20 


1.74 


0.96 


0.95 


C4.5„ 


0.10 


0.05 


15 


1.16 


1.58 


0.99 


0.97 


CV5A 


0.12 


0.02* 


90 


1.55 


1.50 


0.91 


0.75 


CP5.2* 


0.10 


0.01 


590 


1.33 


2.20 


0.83 


0.76 


U6A 


0.33 


0.29* 


8 


1.10 


1.60* 


1.00 


1.00 


U6.2 


0.20 


0.01* 


6 


1.10 


1.56* 


1.00 


1.00 


U6.3 


0.06 


0.14 


9 


1.11 


1.50* 


1.00 


1.00 


U6A 


0.00 


0.09 


8 


1.11 


1.60* 


1.00 


1.00 


54.2 


0.00 


0.09 


506 


2.13 


0.29 


0.99 


0.98 


54.3 


0.10 


0.26 


50 


1.50 


0.97 


0.98 


0.98 



5. The products of collisionless collapse 

Table|2]lists for each simulation the following information: the 
relative mass loss for the end-products AM = (Mo - M)/Mq, 
the relative conservation of the global quantity Q, with AQ = 
I Qo ~ Q\/Qo referred to v = 1/2 unless otherwise noted, the 



concentration C p = p(0)/p(tm) of the end-products in terms 
of the ratio of the central density to the density at the half- 
mass radius, the global anisotropy parameter k = 2K,-/Kt, the 
anisotropy radius (defined by the relation a(r a ) = 1) relative 
to the half-mass radius r a /rM, and the final shape parameters e 
and 77. All quantities are referred to the final system of bound 
particles. 
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Fig. 4. Correlations between final concentration C p and initial 
virial ratio u (left), and between final ellipticity 77 and final 
global anisotropy 2K, /Kt (right). Symbols mark the various 
sets of simulations as follows: filled squares for CI, filled tri- 
angles for C2, open pentagons for C3, open squares for C4, 
crosses for U, and stars for S . 



5.1. General properties 

From the results reported in Table |3 we may infer some em- 
pirical trends. In particular, here we focus on: (1) central con- 
centration; (2) anisotropy content; (3) deviations from spher- 
ical symmetry; (4) mass loss. Density and anisotropy profiles 
will be discussed and compared with our theoretical models in 
Sectional 

We have run two series of simulations (type CI and C2) 
for which the initial particle positions and velocities are kept 
fixed within each series, except for a constant scaling factor in 
the velocities able to lead to different values of u (from 0.05 
to 0.275 for CI and from 0.06 to 0.24 for C2). This procedure 
thus allows us to explore the role of the initial virial ratio by 
keeping all other conditions strictly fixed. 

The central concentration resultin g from the co llapse is ex- 
pected to correlate with u. lLondrillo et alJ dl99lh proposed a 
simple criterion to set an upper limit to the expected value of 
the central concentration by imposing the conservation of the 
maximum density in phase space. They argued that, for the col- 
lapse of an initially homogeneous system, the central concen- 
tration, measured in terms of the ratio tm/?"o.i (of the half-mass 
radius to the radius of the sphere containing one tenth of the 
total mass) should scale as l/u. Our CI and C2 simulations fol- 
low qualitatively the proposed trend. However, since relaxation 
is incomplete, it is natural to find that other factors, in addition 
to the value of u, can contribute to determine the properties of 
the final states. In fact, if we do not restrict our attention to the 
CI and C2 sequences only and consider instead the entire set 
of simulations, we see that the correlation between u and C p 
becomes weaker (see Fig.©. 

Differently from CI and C2, the sets of C3 and C4 simu- 
lations, starting from different spatial configurations (different 
number and size of clumps, different seed in the random num- 
ber generator), allow us to study other possible correlations, 
in particular those between initial and final concentration and 
between final concentration and initial deviations from spher- 
ical sy mmetry; the latter correlation was noted bv lBoilv et alJ 
(2002), starting from homogeneous spheroids. Again, if we in- 



clude the entire set of simulations, the correlations that we find 
are, in general, relatively weak. 

The final global anisotropy of the simulations (see the quan- 
tity k in Table |3 is also weakly correlated with u, with larger 
values of u preferentially associated with lower levels of radial 
anisotropy. The series CI has a systematic, but curiously non- 
monotonic trend, while C3 and C4 show that other factors, in 
addition to u, are important. 

As to the shapes of the products of collisionless collapse, 
we note a relatively strong correlation (see Fig.|4j between the 
final shape (as measured by rf) and the final level of global 
anisotropy (as measured by 2K r /Kj). This is likely to be re- 
lated to the action of the radial orbit instability during collapse. 
In particular, for the C2 series lower values of u lead to more 
anisotropic and more flattened end products; the effect in the 
CI series is less pronounced. Of course, the issue of the final 
shapes produced by collapse has been addressed by several in- 
vestigations in the past, especially with the hope to establish 
whether related dynamical mechanisms can account for the ob- 
served morphologies of elliptica l galaxies (for simulatio ns in 
the co smological context, see lWarren et al J 19921 see also lUdrvl 
119931) . 

Initial conditions with a small number of clumps, as con- 
sidered in our paper, often show significant deviations from 
spherical symmetry (from Table^we see that 770 can be as low 
as 0.7). Curiously, the final value of rj may even slightly exceed 
the value of 770, thus showing that collapse may sometimes push 
the system toward spherical symmetry, not necessarily away 
from it. 

Collisionless collapse can produce significant amounts of 
unbound particles and consequently give rise to mass loss. This 
effect is particularly severe in the cases where th e collapse orig- 
inates from a homogeneous sphere (see also lLondrillo et alJ 
Il99lh : here the system may lose up to one third of the mass 
(see run t/6.1). dumpiness appears to have a stabilizing ef- 
fect with respect to mass loss; in fact, the mass lost is less than 
7% even for run Cl.l characterized by u — 0.06. On the other 
hand, symmetrized clumpy initial states, of type S, are also 
found to evolve with limited mass loss. [Since the nature of the 
gravitational forces is mainly radial for both the collapsing ho- 
mogeneous spheres (U simulations) and symmetrized clumpy 
configurations (S ), the different amounts of mass loss might be 
related to the different radial density distributions for the two 
types of run. In fact, the effect of superimposing several clumps 
of particles creates a density profile decreasing approximately 
linearly with radius.] 

5.2. The role of the radial orbit instability 

Spherical stellar systems with an excess of radial orbits 
(2K,/Kt > 1.7 + 0.25) are expected to be unstable and to 
evolve rapidly, on the dynamical time-scale, into ellipsoids; 
the precise value for the onset of the radial-orbit instability 
depe nds on the detailed structur e of the system co nsidered 
JPolvachenko & Shukhmar]ll98lt see |PalmeJfl99l and ref- 
erences therein). The radial orbit instability is thought to act 
efficiently during collisionless collapse and is then argued to 
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be the leading mechanism that makes cold and spherical ini- 
tial configurations evolve into generally triaxi al configurations 
jAguilar & Merritllll99fi |Polvachenkolll992h . The instability 
may also be responsible for a reduction of th e value of the cen- 
tral c oncentration reached during collapse jMerritt & AguilaJ 
Il985h : in fact, the evolution of concentrated anisotropic sys- 
tems into ellipsoids is accompanie d by a drastic softening of the 
density distribution JStiavelli & Sparkelll99ll) . As is the case 
for many other unstable systems, evolution tends to remove the 
source of instability and thus, in our case, to decrease the ini- 
tial excess of radial orbits. Therefore, the threshold of instabil- 
ity should provide an upper limit to the global anisotropy of 
objects produced by collisionless collapse. 

Our simulations largely confirm the general valid- 
ity of this picture and the g eneral applicability of the 
IPolvachenko & ShukhmarJ dl98lh criterion. In particular, sim- 
ulations C2.3 and C2.4 are characterized by a value of k > 1 .7 
and lead to more flattened configurations than C2.1 and C2.2. 
Also the drop in the central concentration in simulation CI. 10 
with respect to C 1 .9 might be related to the action of the radial- 
orbit instability. Most of the end states are characterized by rel- 
atively high anisotropy (generally k > 1 .5 and values around 
1 .7 are not infrequent) and thus it seems that evolution tends to 
prefer a state very close to the stability boundary (as studied for 
the / (v) family of models in Paper I by means of an extensive 
set of simulations). [An interesting finding is that symmetrized 
initial conditions, although artificial, can lead to spherical fi- 
nal states still able to sustain a large number of radial orbits 
(k =s 2.1 for simulation 54.2)]. 

5.3. Angular momentum mixing 

Simulations with homogeneous initial conditions generate 
quasi-equilibrium final configurations that not only suffer from 
significant mass loss, but also exhibit unusual features in their 
anisotropy profiles (see Sect. l7.3l and Fig. [HJ. 

If the degree of symmetry in the initial conditions is ex- 
cessive, little room is left for relaxation in the (E, J 2 ) phase 
space even if the process itself may be violent and lead to mass 
shedding. This is confirmed by the fact that little or no mixing 
is observed in the single-particle angular momentum distribu- 
tion for homogeneous s imulations, as reported in Fig. |5] (see 
also lMav & van Albadal 19841) . In fact, if the system evolves re- 
maining close to spherical symmetry, the conservation of single 
particle angular momentum imposes severe constraints on the 
dynamical properties of the end-state of the collapse. On the 
other hand, a certain degree of dumpiness, even if limited to 
either position or velocity space, leads to angular momentum 
mixing. This is confirmed by two test simulations, CV5.1 and 
CPS. 2*, where mixing indeed turns out to be quite efficient and 
leads to J relaxation much like in the left panel of Fig. [5] (see 
also Appendix). 

Clumps thus help the system reach a "universal" final state 
from a variety of initial conditions, which can explain the sim- 
ilarity of the density profiles observed in the final products of 
collapse simulations (see Sect.fjji. 




1 2 3 1 2 3 



J(0) J(0) 

Fig. 5. Scatter plot (final vs. initial values) for the single- 
particle specific angular momentum. Comparison between a 
clumpy simulation (run C4.2; left panel) and its symmetrized 
version (run 54.2; right panel). Units for J are pc 2 /yr, see 
Sect. E| 

5.4. Dependence on the degree of dumpiness 

A few simulations with a large number of clumps (400 in C4.5 
and 80 in C4.3) and a spatial filling factor above unity confirm 
that, in the limit of large Nc, the evolution of the system ap- 
proaches that of collapse simulations based on homogeneous 
conditions, with end-states characterized by a flat core and a 
low anisotropy content. A number of clumps of order 10 to 20 
thus seems to be optimal for an efficient violent relaxation. 

Even when limited to either position or velocity space, 
dumpiness can be important and still lead to end-states with 
general properties similar to those of the standard clumpy sim- 
ulations considered in this paper (see CV5.1 and CPS. 2 entries 
in TablesGEland Sect. 1731. 

We also studied the dependence of the results of collision- 
less collapse on the spatial filling factor of the clumps. To 
do this, we took advantage of the ability of GyrFalcON to 
deal with systems with different scales and ran a simulation 
(C4.4) initialized with 80 small cold clumps (i.e. with a radius 
Rc = 2.8 kpc distributed in a sphere of radius 40 kpc). For 
this simulation, evolution basically occurs in two stages, with a 
first collapse in which strongly bound structures are formed in 
a very short time, followed by subsequent merging (see Fig.|6}. 
Interestingly, the outcome of this simulation is highly isotropic 
(a » out to the half-mass radius) and very concentrated. We 
will see (Sect. 17. 5> that, even in this case, the density profile 
remains very well represented by / (v) models (and by the R 1 ^ 4 
law). For C4.4, after several tens of dynamical times, there re- 
main traces (remnants) of the more strongly bound clumps, or- 
biting within the smooth system. 

6. Conservation of Q 

We recall that the quantity Q (for a discrete system of N par- 
ticles Q — 2;=i,;v(-A7|£/l 3 ^ 4 ) v ; see aiso Sect. |3 has been intro- 
duced for the description of conditions in which partial violent 
relaxation occurs, where it is argued that information about the 
initial state is basically lost, except for an approximate conser- 
vation of a third quantity (in addition to total energy and total 
number of particles). Therefore, it would be wrong to invert 
the argument and imagine that, by itself, the conservation of a 
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Fig. 6. Evolution of the virial ratio during run C4.4, character- 
ized by the presence of many small cold clumps. The insert 
box zooms in on the evolution at the beginning of the simula- 
tion, when a first collapse occurs, followed by an expansion of 
the clumps while collapsing toward the center of mass of the 
system. The virial ratio in the final stages of the simulation is 
slightly above unity because of mass loss. 



quantity such as Q is equivalent to the picture of incomplete 
violent relaxation. In particular, we note that, judging from our 
set of simulations, Q is well conserved for homogeneous initial 
conditions, both in the velocity and position space. However, 
this is less relevant to our goals, since homogeneous conditions 
do not allow for mixing and violent relaxation at the level of 
angular momentum space to operate properly. Therefore, it is 
not surprising to find that the end-products of simulations with 
homogeneous initial conditions tend to be less well represented 
by the / (v) models, in spite of their relatively good conservation 
of g. 

In this Section we will show that the issues involved in 
the conjectured conservation of Q and the indications obtained 
from our simulations are complex. Therefore, it would be 
pointless to continue further in this direction, looking for a bet- 
ter definition of what might be defined as "acceptable degree 
of conservation" or searching for other quantities that might 
be conserved better than Q. Instead, to make a decisive test 
about the merits of our approach, we should take the models 
that have been constructed (by means of the Ansatz of the Q- 
conservation) and compare them in detail with the results of 
collisionless collapse obtained from our simulations. Such test 
will be addressed in the following Sect.0 



6.1. The "observed" conservation 

The value of Q, computed with v - 1/2, is approximately con- 
served for a wide range of initial configurations. By approxi- 
mate conservation we mean that AQ < 0.5, although in some 
cases we have conservation as good as AQ « 0.01. As a gen- 
eral rule, Q is better conserved if the initial virial ratio is not 
too low. 

A curious property is that all clumpy simulations appear 
to lead to the same value of Q, with a scatter on the order of 
10% (see Table[3j. [I.e. the scatter is less than the mean devi- 
ation from exact conservation, around 20 - 30%.] This result 
can be interpreted, at the level of the simulations, by consid- 
ering that, independently of the specific details of the initial 
clumpy conditions, the large scale structure of the end products 
of the simulations is very similar, with respect both to physi- 
cal scales (constrained by the conservation of mass and energy 
in the collapse) as well as dimensionless dynamical properties 
at large radii (see Sect. 7). In addition, the fact that the values 
of (M, E to t, Q) realized at the end of the simulations are ap- 
proximately constant is consistent with the fact that the best-fit 
models do not exhibit wide variations in the values of *F and 
v (cf. Table 3 and the discussion of parameter space given by 
Bertin & Trenti 2003, Sect. 3). 

Strict conservation is not meaningful, for a number of rea- 
sons. Indeed, during collisionless collapse even the total num- 
ber of particles N and the total energy E m are not conserved, 
if we refer these quantities to the fin al set of bound particles; it 
was noted ( Stiavell i~& Bertinll 19871) that the non-conservation 
of Q actually correlates with the non-conservation of N and 
E tot . A simple argument also warns us that the conservation of 
Q should not be meant to apply to all conditions. The reason 
is that, if we refer to the proposed definition, Q cannot be con- 
served in the limit of an infinitely cold collapse. In fact, for an 
infinitely cold collapse (i.e. for u — > 0, with the stars kept at 
fixed initial positions), at the beginning of the simulation we 
would have Q — > (because the single-particle angular mo- 
menta vanish, in the limit of vanishing initial velocities, while 
the single-particle binding energies remain at a finite value). On 
the other hand, at the end of the simulation, the formation of a 
quasi-isotropic core with finite kinetic energy content requires 
that the final value of Q be finite. 

Furthermore, the quantity Q is referred to an ideal case 
characterized by spherical symmetry, while, as noted earlier, 
both the initial and the final configurations in our simulations 
of collisionless collapse can exhibit significant deviations from 
spherical symmetry. To get an estimate of changes of Q associ- 
ated with deviations from spherical symmetry, we have consid- 
ered a (1/2; 3) / (v) model, unstable against the radial orbit in- 
stability (Paper I), and let it evolve; the final quasi-equilibrium 
state is characterized by e * rj « 0.73 and is associated with 
a change AQ = 0.12. Similar changes are observed by stretch- 
ing artificially an / (v) model to a non-spherical geometry, with 
e = 1 and rj « 0.7. But these changes are given for comparison 
only, since they are not related to conditions in which violent 
relaxation takes place. 
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Fig. 7. Conservation of the general functional Q (see Eq. m 
a typical simulation (C4.1). 

6.2. General polynomial dependence 

Although aware of the fact that we should not really look for 
quantities conserved exactly during collisionless collapse, we 
decided to test the paradigm of Q conservation further by con- 
sidering the more general class of Q functionals, defined, for a 
system of N points, as: 



2 = Z 



~t \Ei\ 



(1) 



where v\ and V2 are free parameters. We explored the parameter 
space — 1 < vi < 1 and -1 < V2 < 1. The functional Q used to 
construct the / (v) models corresponds to the condition v = v\ = 
V2 > 0, which guarantees the desired asymptotic behavior for 
the associated density p ~ r~ 4 at large radii. 

We studied the change in the value of this functional com- 
puted at the beginning and at the end of a typical simulation 
(10 5 particles in 10 cold clumps, run C4.1; see Fig.0. If we 
focus on the v\ — V2 = v condition, the best conservation would 
be attained for low values of v. 

7. Fit with the f (y) models 

We first fit the density and the pressure anisotropy profiles, p(r) 
and a(r), of the end-products of our simulations by means of 
the / (v) family of models. The phase space properties of the 
best-fit model thus identified are then compared with those of 
the end-products of the simulations. 

Smooth, angle -averaged simulation profiles are obtained 
by binning the particles in spherical shells and averaging over 
time, based on a total of 20 snapshots taken from t = 64 to 
t = 80, at an epoch when the system has already settled down 
in a quasi-equilibrium configuration. For the / M models, the 
parameter space explored is that of an equally spaced grid in 
(v, with a subdivision of 1/8 in v, from 3/8 to 1, and of 
0.2 in y ¥, from 0.2 to 14.0 (corresponding to the grid of models 



studied in Paper I). The mass and the half-mass radius of the 
models are fixed by the scales set by the simulations. 

A minimum-^ 2 analysis is then performed, with error bars 
estimated from the variance in the time average process used 
to obtain the smooth simulation profiles. A critical step in this 
fitting procedure is the choice of the relative weights for the 
density and the pressure anisotropy profiles. We adopted equal 
weights for the two terms, checking a posteriori that their con- 
tributions to x 2 are of the same order of magnitude. 

Table 3. Best fit / (v) models for the set of high resolution runs 
(series C2 and C3). The various columns give: run identifier, 
model identifier, mean value of the absolute relative deviations 
from the density of the simulations, mean value of the absolute 
deviations in the pressure anisotropy profile, mean value of the 
absolute relative deviations in the energy distribution, and final 
value of Q. 





f M 


<|Ap/p|> 


(|Aa|> 


(\AE/E\) 


Q 


C2.1 


(l/2;4.8) 


0.11 


0.07 


0.23 


1.24 


C2.2 


(l/2;4.8) 


0.11 


0.06 


0.22 


1.33 


C2.3 


(5/8;5.0) 


0.12 


0.06 


0.21 


1.35 


C2.4 


(7/8,5.6) 


0.14 


0.08 


0.23 


1.33 


C3.1 


(3/8 ;5. 6) 


0.10 


0.22 


0.18 


1.33 


C3.2 


(3/8;5.4) 


0.11 


0.19 


0.22 


1.26 


C3.3 


(l/2;5.2) 


0.17 


0.16 


0.20 


1.64 


C3.4 


(5/8 ;5. 4) 


0.12 


0.05 


0.18 


1.40 


C3.5 


(l/2;6.2) 


0.09 


0.20 


0.15 


1.35 


C3.6 


(3/8;5.2) 


0.13 


0.05 


0.20 


1.35 



7.1. Density profiles 

Since the half-mass radius tm and the total mass M are kept 
fixed in the fitting procedure, we are left with two degrees of 
freedom (i.e., the dimensionless parameters v and *F). In prac- 
tice, given the general behavior of the density profile of the 
models (see Fig. 3 in Paper I), at large radii the freedom 
in the fit is limited. Therefore, the excellent match at large 
radii to the density profile of the end-products of the simula- 
tions demonstrates that the / (v) family has been constructed 
on solid physical grounds. Different values of (v, corre- 
spond to different shapes of the inner potential well and of 
the anisotropy profile. As exemplified by Figs. I9I1 II the den- 
sity of the final systems produced by the high resolution set 
of simulations (C2 and C3) is well represented by the best- 
fit f v) profile over the entire radial range, from 0.1 to 10 half 
mass radii. The fit is satisfactory not only in the outer parts, 
where the density falls by nine orders of magnitude with re- 
spect to the central regions, but also in the inner regions. The 
mean absolute relative deviation between simulations and mod- 
els «|Ap/p|> = (1/Ay I^Ji \Psim(n) -p m0 dei(n)\/psim(ri)), com- 
puted over this extended radial range, is usually around 10% 
(see Table here N g represents the number of radial grid 
points. 



M. Trenti et al.: A family of models of partially relaxed stellar systems 



11 



With a similar procedure, we have studied the end-products 
of simulations characterized by different numbers of particles 
and clumps (CI and C4). No significant changes in the quality 
of the fits are found if we focus on simulations characterized by 
clumpy initial conditions (with the possible exception of those 
run with N c > 80) . 

7.2. Projected density profiles 

The end-products of collisionless collapse are known to be 
characterized by projected density prof iles generally well fitted 
by the R 1 ^ 4 law Jde Vaucouleursll 19481) . provided that the col- 
laps e factor is large (i . e., that the initial vir ial ratio u is small; 
see Ivan Albadalll982l lLondrillo et alJll99lh . With our set of 
simulations we confirm this result and we extend it by means 
of the / (v) models. 

The successful comparison between models and simula- 
tions is interesting because, depending on the value of u, some 
simulations lead to configurations that exhibit deviations from 
the R l l 4 law. In these cases, the density profile projected along 
the line-of-sight is characterized by an R ^" behavior with 
n + 4. For example the C2.4 simulation, which starts with a 
low collapse factor, has a best fit index n m 3, while the the 
simulation C3.1, which has a large collapse factor, is best rep- 
resented by a profile with n « 5. Yet, these systems turn out to 
be all well fitted by the / (v) models. Therefore, the family of 
models that we have identified might also be useful to describe 
systematic structural changes in galaxies, in the f ramework of 
the p roposed weak homology of elliptical galaxies ( B ertin et aD 
l2002l) . 

7.3. Pressure anisotropy profiles 

In our simulations the pressure anisotropy profiles follow the 
general trend expected for the process of collisionless collapse. 
In particular, the final configurations are characterized by an 
isotropic core, with a ~ 0, while the outer regions have a 
strongly radially biased anisotropy (up to a — 2). The transi- 
tion region (a ~ 1) is located around the half-mass radius (see 
column r a /rM in Tabled Higher values of 2K r /Kr are associ- 
ated with lower values of r a /rM- For clumpy initial conditions 
(with the possible exception of those run with Nc > 80), the 
anisotropy profile a(r) is a monotonic increasing function of 
the radius. A curious feature is found for the results of collapse 
of uniform spheres (runs U). Here (see Fig.|S|l the core is ba- 
sically isotropic, with the region around the half-mass radius 
exhibiting an excess of tangential orbits (up to a ~ -0.4). In 
the outer parts, but with a very sharp transition, the pressure 
profile becomes radially biased. In correspondence to the dip 
in a, where a < 0, we note a clear feature in the density profile 
(see Fig.|8|l. Uniform spheres initialized with a very small par- 
ticle number (iV < 10 4 ) do not show this behavior; for them the 
pressure anisotropy rises quite regularly, although the profile is 
significantly affected by Poisson noise. 

In conclusion, for all the clumpy C runs (again, we 
should mention, with the possible exception of those runs with 
Nc > 80), the anisotropy profile is represented extremely 




-1 

Log (r/r.) 




Log (r/r u ) E 



Fig. 8. Density and anisotropy profiles (left frames) and energy 
density distribution (right frame) for simulation U6.2, starting 
from a homogeneous sphere. Note that in the vicinity of the 
half-mass radius the pressure anisotropy is tangentially biased. 

well by our models, with a mean absolute error ((IAcI) = 
(W«)£i=i \a S im(n) - a mode i(n)\) typically around 0.1 but of- 
ten as low as 0.05 (see Table|3}. 

To some extent, the final anisotropy profiles for clumpy ini- 
tial conditions are found to be sensitive to the detailed choice 
of initialization. In other words, runs starting from initial con- 
ditions with the same parameters, but with a different seed in 
the random number generator, give rise to slightly different pro- 
files. In any case, the agreement between the simulation and the 
model profiles remains very good (see Figs. 1911 TV 

7.4. Comparison at the level of phase space 

At the level of phase space, we have performed two types of 
comparison, one involving the energy density distribution N(E) 
and the other based on N(E, J 2 ). The chosen normalization fac- 
tors are such that: 

M = j N(E)dE = j N(E, J 2 )dEdJ 2 . (2) 

The energy distributions N(E) that we find (see Fig. 1911 1\ . 
qualitatively similar to those obtained in earlier investigations 
(see Fig. 2 in Ivan Albadal 19821 and Fig. 10 in lUdrvlll993l) . are 
characterized by an approximate exponential behavior at low 
energies (N(E) oc exp(-aEJ) with a rapid cut-off near the ori- 
gin, which is argued to go as \E\ 5 ^ 2 bec ause the potential is 
Kepleri an in the outer par ts dUdrvll 19931 see also the discus- 
sion bv traffdll987l and bv lBertin & Stiavellilll989l) . The final 
states of the simulations also show the presence of particles 
with positive energy, escaped from the system. 
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Fig. 9. Comparison between the C3.5 simulation and the best- 
fit f iy) model (1/2; 6.2). The top left panel represents the den- 
sity as measured from the simulation (error bars) and the best- 
fit profile (line). The top right panel gives the residuals from 
the fit. At the bottom left, the anisotropy profile of the simu- 
lation (error bars) is compared with the best-fit profile (line); 
the bottom right frame illustrates the energy density distribu- 
tion N(E). The density p and the single-particle energy E are 
given in code units. 



In Fig. [5] we plot the final energy density distribution for 
the simulation run C3.5 with respect to the predictions of the 
best-fit model identified from the study of the density and pres- 
sure anisotropy distributions. Similar plots are given in the fol- 
lowing figures for other simulations. The agreement is very 
good {{\hE\) as 0.2, see Table [3}, especially for the strongly 
bound particles. In particular, this means that we are correctly 
describing the innermost part of the system. The energy dis- 
tribution for less bound particles (i.e. those associated mostly 
with the outer parts of the system) is less regular and sometimes 
presents a double peak (e.g., see Fig.1101. which obviously can- 
not be matched in detail by our models. This is an interesting 
example of the way some memory of the initial state can be 
preserved (the extra-peak is indeed related to the initial distri- 
bution of binding energies) and a direct sign of the incomplete- 
ness of violent relaxation. 

Finally, at the deeper level of N{E, J 2 ), simulations and 
models also agree rather well, as illustrated in the four panels 
of Figs. 1121141 For the cases shown, the distribution contour 
lines are in good agreement in the range from E m - m to E x -4; 
however, the theoretical models show a peak located near the 
origin, not present in the simulations, which is related to the 
Jacobian factor arising from the transformation of the dis- 
tribution function from the (x, w) to the (E, J 2 ) space. 



Fig. 10. Comparison between the C2. 1 simulation and the best- 
fit f v) model (1/2; 4.8), shown as in Fig.|9] 




Fig. 11. Comparison between the C3.4 simulation and the best- 
fit / (v) model (5/8; 5.4), shown as in Fig.|9] 



7.5. An additional test to characterize clumpy initial 
conditions 

As an additional test to characterize the detailed effects of 
dumpiness, we studied the end-products of the CV5.1 and 
CP5.2* simulations, by comparing them with the / (v) models. 

Although these two runs start from initial conditions rather 
different from our standard choice (cf. C1-C3), being homo- 
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Fig. 12. Final phase space density N(E, J 2 ) (left column) for Fig. 14. Final phase space density N(E, J 2 ) (left column) for 
the simulation C3.5, compared with that of the best fitting the simulation C3.4, compared with that of the best fitting 



(1/2; 6.2) f v > model (right column). 



(5/8; 5.4) f v> model (right column). 
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Fig. 13. Final phase space density N{E, J 2 ) (left column) for 
the simulation C2.1, compared with that of the best fitting 
(1/2; 4.8) f {v) model (right column). 

geneous either in position (CV5.1) or in velocity (CP5.2*) 
space, we note that they can be fitted very well by our family 
of models: (3/4; 5.4) for CV5.1 and (1;6.2) for CP5.2* (with 
(\Ap/p\) as 0.1). The good match at the level of the anisotropy 
profile a(r) and of the single-particle energy distribution also 
confirms, as discussed in the Appendix, that the requirement of 
dumpiness in phase space is a well posed characterization of 
the initial conditions. The two runs have the following behavior 
with respect to the 2- conserva tion: AQ = 0.02 and final value 
Q = 1.40 for CV5.1; AQ = 0.01 and final value Q = 1.26 for 
CP5.2*. 

In passing, we note that the C4.4* simulation, character- 
ized by very small clumps, leads to a concentrated final density 
profile that is well reproduced by the (1; 9.2) / (v) model (with 
(\Ap/p\) * 0.15). 



7.6. Separate fits to density and anisotropy profiles by 
means of simple analytic functions 

Simple analytic descriptions of density profiles and, separately, 
of anisotropy profiles are often used in stellar dynamics, with- 
out a specific physical scenario of galaxy formation. For the 
density profile we may refer to: 

(3 - y)M 



4n ryir + rof-y' 



(3) 



where < y < 3 is a free par ameter, a nd M and ro are a mass 
and length scale respectively llDehnenlll993h . As discussed in 
Paper I, it is no surprise to find that the case y = 2 Jjaff el 19831) 
captures the general properties of the density profile obtained 
by the simulations at the 20% level. Curiously, when we fit the 
density distribution of some simulations by means of Eq. l|3}, 
the best fitting index y is very low y w 0.1 (see Fig. I15i . 

Similarly, for the anisotropy profile one might resort to the 
analytic distribution 



o(r) = 2- 



1 ' 
r- 



(4) 



with r a being a free scale llMerritJl 985). As shown in Fig. 1151 
the typical shape of the anisotropy profile reached at the end of 
the simulations is different. 

8. Discussion and conclusions 

In this paper we have concentrated on nearly spherical, one- 
component stellar systems. As is well known, in spite of these 
restrictions, the equations of stellar dynamics allow for almost 
complete freedom in the construction of self-consistent dynam- 
ical models, with the only requirement that they should be sup- 
ported by a positive definite (but otherwise arbitrary) function 
of E and J, as a distribution function in phase space. Therefore, 
the full range of self-consistent one-component spherical stel- 
lar dynamical models is enormous. Most likely, the majority of 
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Fig. 15. Density profile (left), fitted using Eq. Q with y = 
0.11, and anisotropy profile (right), fitted with Eq. ®, for 
the simulation C3.4. Compare to the fit with the (5/8 ;5. 4) / (v) 
model in Fig. II II 

these models have little to do with the systems that have been 
realized in nature. The main idea at the basis of the present 
paper is to combine clues from A^-body simulations and from 
statistical arguments so as to pinpoint, among the enormous 
variety of in principle acceptable dynamical models, those few 
that, because of their physical justification, have a chance of 
matching the properties of interesting classes of numerical sim- 
ulations and of observed stellar systems. 

Some interesting clues had been noted earlier. With the 
aim of summarizing the main properties of incomplete vio- 
lent relax ation during collisionless collapse, it was discovered 
JStiavelli & Bertinlll987l) that, by arguing that a third quantity 
Q (in addition to total energy and number of stars) should be 
included among the relevant constraints in the extremization of 
the Boltzmann entropy, the most probable and thus physically 
justified distribution function / (v) leads to models that are in 
general qualitative correspondence with the products of colli- 
sionless collapse found in numerical simulations and with the 
observed luminosity profiles of bright elliptical galaxies. 

In the present paper we have demonstrated that the / (v) 
models are able to match in surprising quantitative detail the 
results of our numerical simulations. At the same time, the 
models exhibit projected density profiles that are well repre- 
sented by the R 1 ^" law (generally with 11*4; the residuals from 
the fit are within 0.1 magnitudes in a radial range from 0.1 to 
10 effective radii; see also Paper I). Therefore, we have demon- 
strated that the / <v) models, as well as the end products of the 
collapse simulations, are relevant to the description of the stel- 
lar distribution of elliptical galaxies. Such correspondence is 
even more remarkable, if we recall that, from the results estab- 
lished in the last decades, dark matter should play a dominant 
role in the structure of galaxies, while our approach neglects, 
so far, some important ingredients among which the presence 
of a massive, possibly diffuse dark halo. 

Independently of stellar dynamical modeling, our simula- 
tions have shown that clumpy initial configurations allow for an 
efficient re-distribution of the angular momenta of the individ- 
ual particles during collapse: such efficient phase space mixing 
is precisely the main condition required for a successful appli- 
cation of the statistical arguments that lead to the construction 
of the / M family of distribution functions. In the past (e.g., see 



van Albadall982tlMay & van Albadall984lMerritt & AguilaJ 
1985; Londril loetalJll99ll) . it has been noted that cold col- 
lapses, within a wide range of initial density profiles, gener- 
ate quasi-equilibrium systems with approximate R 1 ^ 4 profiles. 
Here we confirm that the best match to approximate R 1 ^ 4 pro- 
files is obtained from initiallly clumpy configurations. It thus 
appears that collapses starting from artificially uniform and 
spherically symmetric initial conditions retain too much mem- 
ory of the initial conditions and are unable to evolve into a uni- 
versal density distribution. Therefore, it is interesting to find 
that precisely those initial conditions that look more plausible 
and realistic from the physical point of view lead to end prod- 
ucts able to match the stellar distribution of observed systems 
in detail. We may then conclude that collisionless collapse from 
clumpy initial conditions followed by violent relaxation is in- 
deed a formation mechanism relevant to elliptical galaxies. 

If we now take the point of view of stellar dynamical mod- 
eling and examine the foundation of the / (v) family of models, 
we note that many collapse simulations show ^-conservation 
at the 20% level or better (e.g., Cl.l, C2.1 and C3.4). But it 
is even more surprising to find that the end products can be 
fitted so well by the / (v) models. Such good fits make it clear 
that the assumption of Q conservation narrows down the very 
wide range of self-consistent dynamical models to precisely 
those few systems whose properties match both observed sys- 
tems and the end products of collisionless collapse. One must 
conclude that the value of the ^-conservation assumption goes 
beyond mere "physical plausibility" and "mathematical conve- 
nience": it does serve as a sound physical basis for the construc- 
tion of dynamical models of partially relaxed stellar systems. 

We should emphasize that such detailed quantitative cor- 
respondence with observed systems and with the end products 
of collisionless collapse comes as a complete surprise, because 
the two parameters that can be varied within the / (v) family of 
models (i.e., v and T) leave very little freedom with respect to 
density and anisotropy profiles (see Paper I). Especially note- 
worthy are not only the match of the density profile over nine 
orders of magnitude but also the excellent agreement of the ve- 
locity anisotropy profiles between the / (v) models and several 
end products of collapse from clumpy initial conditions (see 
Figs. I91TT1 and Tabled. 

Yet, one cannot claim that the / (v) models give a fully sat- 
isfactory description of the phase space structure of systems 
produced via incomplete violent relaxation. In fact, the associ- 
ated N{E, J 2 ) distribution is characterized by singular behav- 
ior near the origin in the (E, J 2 ) plane, which is not present in 
the end-states of the simulations. In spite of such discrepancy 
between models and end-products of the simulations, the inte- 
grated properties (e.g.,N(E), a(r) and p(r)) are very well repro- 
duced. This confirms the fact that a variety of different distribu- 
tions in phase space can lead to the same integrated properties. 
In this respect, it appears that, if we refer to the extreme outer 
parts of the system (with r » r^, and £ — > 0) the previously 
studied f x models ( Be rtin & Stiavellill 19841) . with their regu- 
lar distribution function f(E, J 2 ) « \E\ 3 ^ 2 at low values of \E\, 
might still have an advantage over the / M models. 

Another interesting (although partly known) result of the 
present paper is that the velocity distributions of the end prod- 
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ucts of the collapse simulations and of the best fitting mod- 
els possess, in many cases, a rather strong radial anisotropy. 
In some of the collapse simulations we see clear signs that 
the radial-orbit instability has been active (as indicated by the 
correlation between final ellipticity rj and anisotropy content 
2K r /Kr\ cf. Fig.@J, resulting in end products that are close to 
the threshold for the onset of the radial-orbit instability. In gen- 
eral, systems that are unstable with respect to the radial-orbit 
instability should evolve into marginally stable systems (see 
also the study of the unstable (1; 3.2) / (v) model in Paper I). In 
view of the good correspondence between the results of the for- 
mation processes studied in this paper and important observed 
properties of elliptical galaxies, we may argue that ellipticals 
are also likely to lie close to the threshold of radial-orbit insta- 
bility. This would happen if elliptical galaxies, during their for- 
mation process, indeed went through a collisionless phase char- 
acterized by strong radial motions (such as collapse or head- 
on mergers). We plan to better quantify this connection by ex- 
tending the study to two-component models and collapses, also 
starting from a power spectrum of perturbations representative 
of cosmological initial conditions. 

The last remark brings us naturally to one final comment. 
We recall that, since collisionless dynamics is scale-free, the 
results obtained here can also be interpreted as relevant to the 
description of the collapse of dark matter halos. Clearly, since 
we do not include the effects related to the general Hubble ex- 
pansion and we do not initialize our clumpy conditions in terms 
of the power spectrum of perturbations appropriate for a given 
cosmological epoch, a direct comparison between our set of 
numerical experiments and the profiles of dark ma tter halos ob- 
tained in ACDM simulations ( Navarro et alJl997tlMoore et alJ 
1998) would not be justified. Still, our experiments can be con- 
sidered as one example of final equilibrium realizations of a 
dark halo, when initial conditions are varied outside the pre- 
scriptions consistent with the curre ntly accepted cosmological 
framework (see also iLemsonll 19951) . If we now go back to our 
interpretation in terms of the / <v) models, it is noteworthy to 
point out that, although the density profile of the / (v) models 
falls off as 1/r 4 at large radii, in the inner parts that might 
correspo nd to the regions in side the virial radius (for a defi- 
nition see lNavarro et alll997h . the density goes approximately 
as 1/r 3 2 (see Sect. 3.1 in Paper I), that is very c lose to the re- 
ported 1/r 3 value for co smological simulations dNavarro et alJ 
ll997HMoore et all 1 9981) . Since the outskirts of dark matter ha- 
los are "still collapsing", and thus their dynamical conditions 
are different from those under which we derived the / (v) mod- 
els, this agreement appears surprisingly good and suggests fur- 
ther investigations. 
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Appendix A: A quantitative measure of 
dumpiness 

In order to characterize the degree of dumpiness present in the 
initial conditions of our simulations, we may consider, in the 
6-dimensional phase space, the ratio cl = ipf 1 J)l{p^) of the 
mean local density around particles to the mean density. 

We estimate the mean 6-dimensional density in phase space 
(p (6) ) by dividing the number of particles N by the typical to- 
tal volume occupied. Since the large-scale structure in phase 
space is that of a sphere both in position and velocity space 
separately, we compute the total volume as the product of these 
two volumes. Each volume is calculated by assuming that the 
radius of each sphere is equal to the mean distance between two 
randomly chosen particles in the relevant space (position and 
velocity respectively); for example, for a homogeneous density 
distribution inside a sphere of unit radius, the radius determined 
from the adopted procedure would be « 1.03. 

The local density pfl al (required for calculating the aver- 
age used in the definition of cl) is computed by considering 
one particle and by counting the number of neighboring par- 
ticles Niocai within a six-dimensional small sphere of fixed ra- 
dius r s (and thus by assuming an equally weighted norm in the 
phase space for positions and velocities). The scale r s is cho- 
sen in such a way that, on average, a small fixed fraction of 
the total number of particles is enclosed. We set this fraction to 
be £ = {N local) IN ~ 1/250. This choice ensures that we have, 
on average, a high filling factor within the small sphere, so that 
the effects of biases in the local density estimation, arising from 
the coincidence of the center of the local sphere with the coor- 
dinates of a particle, are unimportant (for a discussion on the 
construct ion of unbiased es timators for the local density, see 
also lCasertano & Hull985T) . 

The adopted scale r s also acts as a cut-off scale to the 
dumpiness estimator cl, which is obviously insensitive to fluc- 
tuations at scales smaller than r s . The dependence of the 
dumpiness estimator on £ is illustrated in Fig. IA.ll Eventually, 
diagnostic tools such as cZ(^), as a measure of the initial spec- 
trum of inhomogeneities in phase space, will help us establish 
a bridge toward initial conditions representative of the cosmo- 
logical context (see also comments at the end of Sect. 14. U . 

For our homogeneous initial conditions (simulations of 
type U) the value of the dumpiness estimator is 0.65 < cl < 1, 
depending on the scale considered (cl = 0.72 for £ = 1/250). 
Note that the value of cl can fall below unity, because of bound- 
ary effects. In contrast, for the cold clumpy initial conditions of 
type Cl, C2, and C3 (with 10 and 20 clumps, and spatial filling 
factor N c x R 3 C /R 3 ~ 1.25), at f = 1/250 cl takes on values 
above 30, with typical values around 50 and peaks up to 100. 
For simulation C4.4 (with "small" clumps, and spatial filling 
factor N c x R 3 C /R 3 = 0.027), cl increases to 300. Conversely, cl 
decreases by increasing the number of clumps (down to cl = 15 
for simulation C4.3 with 80 clumps and to cl * 4.5 for simula- 
tion C4.5 with 400 clumps). 

With the numbers quoted above, we see that, at fixed num- 
ber of particles, the dumpiness estimator cl varies with the 
number of clumps Nc used. 



16 



M. Trenti et al.: A family of models of partially relaxed stellar systems 




Fig.A.l. dumpiness estimator cl as a function of £ = 
(Nhcai)/N, for the initial conditions of simulations C4.1 (10 
clumps), C4.3 (80 clumps), and C4.5 (400 clumps). The spa- 
tial filling factor is kept approximately constant (NcXR^/R 3 = 
1.1 - 1.3). The arrow indicates the scale £ = 1/250 at which we 
refer most of our estimates. 

A .1. dumpiness and mixing 

As already anticipated in Sect. I5.3I and I5.4I for an efficient 
angular momentum mixing it is sufficient that dumpiness be 
present either in position or velocity space. In fact, a simulation 
starting from uniform conditions in terms of positions but with 
clumpy structure in velocity space is bound to develop, after a 
few dynamical times, a significant dumpiness in position space 
(see Fig. IA.2L so that the single-particle angular momenta are 
well mixed at the end of the simulation (much like in the left 
panel of Fig. [5}- This result confirms that our choice for quanti- 
fying the dumpiness of a given configuration by looking at the 
six-dimensional phase space is indeed reasonable. 
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